Generated from rmarkdown file: “code/behav/_behav.rmd“. This rmd file executes ‘child’ R-scripts within this directory. These child scripts contain analysis code. Their output is captured and ‘knitted’ back into this report. To run the code, you can either execute this .rmd file (i.e., ‘knit’ to html), or source each child script individually. (These child scripts should be configured to run on their own when sourced directly.)
plot(behav$rt)
abline(h = 250)
abline(h = 3000) ## marks beginning of subsequent trialbehav$is.implausible.rt <- behav$rt > 3000 | behav$rt < 250
grid.arrange(
behav %>%
filter(acc == 1, !is.implausible.rt) %>%
ggplot(aes(rt)) +
geom_density(fill = "slateblue", alpha = 0.3) +
labs(title = "all correct trials") +
theme(legend.position = "none"),
behav %>%
filter(acc == 1, !is.implausible.rt) %>%
ggplot(aes(rt)) +
geom_density(aes(fill = trial.type), alpha = 0.3) +
labs(title = "by congruency") +
scale_fill_brewer() +
scale_color_brewer() +
theme(legend.position = c(0.5, 0.5)),
ncol = 2
)behav %>%
filter(acc == 1, !is.implausible.rt) %>%
full_join(group_by(., subj) %>% summarize(r2norm = qqr2(rt)), by = "subj") %>%
ggplot(aes(sample = rt)) +
stat_qq(alpha = 0.8, size = 1) +
stat_qq_line(size = 0.25) +
facet_wrap(vars(subj)) +
geom_text(
aes(
x = -1.25, y = 2000,
label = paste0("r^2 = ", round(r2norm, 3))
), size = 3, color = "grey50"
) +
theme_minimal(base_size = 10) +
theme(axis.title = element_blank()) +
labs(
title = "QQ rt versus normal",
subtitle = paste0("overall r^2 with normal = ", round(qqr2(behav$rt), 3))
)849971 and 161832 have odd patterns in RT distribution. Strong clustering at 500 ms. This is highly consistent with an artifact we found in other RT data from the same microphone/preprocessing method.
behav %>%
filter(acc == 1, !is.implausible.rt) %>%
ggplot(aes(x = rt, group = subj)) +
geom_density(data = . %>% filter(!subj %in% c("849971", "161832")), color = rgb(0, 0, 0, 0.15), size = 2) +
geom_density(data = . %>% filter(subj %in% c("849971", "161832")), color = "firebrick1", size = 2) +
labs(
title = "subjects with bimodal distributions (likely artifactual) in red"
)behav$is.artifactual.rt <- FALSE
behav$is.artifactual.rt[behav$subj %in% c("849971", "161832") & behav$rt < 500] <- TRUE
behav %>%
filter(acc == 1, !is.implausible.rt, subj %in% c("849971", "161832")) %>%
full_join(group_by(., subj) %>% summarize(r2norm = qqr2(rt)), by = "subj") %>%
ggplot(aes(sample = rt)) +
stat_qq(alpha = 0.8, size = 0.4) +
stat_qq_line(size = 0.25) +
geom_hline(yintercept = 500) +
facet_wrap(vars(subj)) +
theme_minimal(base_size = 10)## make sure to exclude validation set!
behav.rt.aset <- behav %>% filter(acc == 1, !is.implausible.rt, !is.artifactual.rt, is.analysis.group)
fit1.het <- lme(
rt ~ trial.type,
random = ~ trial.type | subj,
data = behav.rt.aset,
weights = varIdent(form = ~ 1 | subj),
control = cl1,
method = "REML"
)
## define "far" outliers as points with resids more extreme than 3 IQR
behav.rt.aset$resid.p <- resid(fit1.het, type = "p")
behav.rt.aset$is.far.out <- farout(behav.rt.aset$resid.p)
## exclude far outliers and refit model
fit1.het.trim <- update(fit1.het, data = behav.rt.aset %>% filter(!is.far.out))
fit1.het.trim.ml <- update(fit1.het.trim, method = "ML") ## for model comparisonsfit1.hom.trim.ml <- lme(
rt ~ trial.type,
random = ~ trial.type | subj,
data = behav.rt.aset %>% filter(!is.far.out),
control = cl1,
method = "ML"
)
(rt.hom.v.het <- anova(fit1.hom.trim.ml, fit1.het.trim.ml))## Model df AIC BIC logLik Test L.Ratio p-value
## fit1.hom.trim.ml 1 6 132710.8 132754.1 -66349.38
## fit1.het.trim.ml 2 54 128508.2 128898.3 -64200.10 1 vs 2 4298.551 <.0001
fit0.het.trim.ml <- update(fit1.het.trim.ml, random = ~ 1 | subj)
(rt.stroopvar <- anova(fit0.het.trim.ml, fit1.het.trim.ml))## Model df AIC BIC logLik Test L.Ratio p-value
## fit0.het.trim.ml 1 52 128591.7 128967.4 -64243.87
## fit1.het.trim.ml 2 54 128508.2 128898.3 -64200.10 1 vs 2 87.53103 <.0001
fit1.het.run.trim <- lme(
rt ~ trial.type * run,
random = ~ trial.type * run | subj,
data = behav.rt.aset %>% filter(!is.far.out),
weights = varIdent(form = ~ 1 | subj),
control = cl1,
method = "REML"
)
summary(fit1.het.run.trim)## Linear mixed-effects model fit by REML
## Data: behav.rt.aset %>% filter(!is.far.out)
## AIC BIC logLik
## 128353.5 128808.6 -64113.76
##
## Random effects:
## Formula: ~trial.type * run | subj
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 117.42723 (Intr) trl.ty run
## trial.typeincon 33.54431 0.006
## run 42.04609 0.315 -0.006
## trial.typeincon:run 11.10237 -0.428 0.079 0.707
## Residual 82.02221
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | subj
## Parameter estimates:
## 107321 123117 130114 132017 138837 141422
## 1.0000000 2.2042702 1.1374693 1.9329246 1.8295601 1.7368040
## 150423 158136 160830 161832 165032 171330
## 2.0389919 1.0513668 4.5850090 5.1845818 2.6912996 1.8198288
## 178243 178950 182436 197449 203418 204319
## 2.8557941 1.3055261 3.2531314 2.1067967 1.1163181 2.5113983
## 300618 317332 346945 352738 448347 562345
## 1.0329181 1.1376613 1.4761585 2.5163872 1.6665109 1.2050780
## 580650 594156 601127 765864 814649 843151
## 3.3349946 1.0921162 1.9364733 1.0036115 1.7891117 2.1310863
## 849971 873968 877168 DMCC1971064 DMCC2442951 DMCC3062542
## 2.7518139 1.2534452 2.1358899 1.1221395 0.7311796 1.7021958
## DMCC5009144 DMCC6418065 DMCC6484785 DMCC6661074 DMCC6671683 DMCC6721369
## 0.8460860 0.9366685 1.7051492 1.3783696 0.9051752 0.9766951
## DMCC7297690 DMCC7921988 DMCC8033964 DMCC8050964 DMCC9441378 DMCC9478705
## 1.3092506 2.0402631 0.8997894 1.6099109 2.4952757 1.5766460
## DMCC9953810
## 1.5094176
## Fixed effects: rt ~ trial.type * run
## Value Std.Error DF t-value p-value
## (Intercept) 787.5251 18.070759 10089 43.58008 0.0000
## trial.typeincon 75.4967 9.060218 10089 8.33277 0.0000
## run 22.3331 7.295533 10089 3.06120 0.0022
## trial.typeincon:run 2.9932 5.017546 10089 0.59654 0.5508
## Correlation:
## (Intr) trl.ty run
## trial.typeincon -0.227
## run 0.047 0.319
## trial.typeincon:run 0.104 -0.725 -0.200
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.8308652 -0.6703528 -0.1488783 0.4980889 4.6330484
##
## Number of Observations: 10141
## Number of Groups: 49
## get unconditional / marginal covariances
m <- rbind(c(0, 1, 0, 0), c(0, 1, 0, 1)) ## contrast matrix
tau.trim <- getVarCov(fit1.het.run.trim)
(tau.trim <- m %*% tau.trim %*% t(m)) ## covariance matrix## [,1] [,2]
## [1,] 1125.221 1154.655
## [2,] 1154.655 1307.351
(r.marginal.trim <- cov2cor(tau.trim)[1, 2]) ## correlation## [1] 0.9520003
svd(getVarCov(fit1.het.run.trim), nv = 0L)$d ## covmat not degenerate (eigvals > 0)## [1] 14005.884684 1670.837724 1126.077069 2.712698
behav <- behav %>% mutate(error = 1 - acc)
fit.error0 <- glmer(
error ~ trial.type + (1 | subj),
behav %>% filter(response.final != "unintelligible", is.analysis.group),
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)
fit.error1 <- glmer(
error ~ trial.type + (trial.type | subj),
behav %>% filter(response.final != "unintelligible", is.analysis.group),
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.093792 (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
(er.stroopvar <- anova(fit.error0, fit.error1))## Data: behav %>% filter(response.final != "unintelligible", is.analysis.group)
## Models:
## fit.error0: error ~ trial.type + (1 | subj)
## fit.error1: error ~ trial.type + (trial.type | subj)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit.error0 3 1729.8 1751.6 -861.89 1723.8
## fit.error1 5 1733.1 1769.4 -861.55 1723.1 0.6921 2 0.7075
fit.error1.run <- glmer(
error ~ trial.type * run + (trial.type * run | subj),
behav %>% filter(response.final != "unintelligible", is.analysis.group),
family = binomial,
control = glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 1E9))
)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.403595 (tol = 0.001, component 1)
summary(fit.error1.run)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: error ~ trial.type * run + (trial.type * run | subj)
## Data:
## behav %>% filter(response.final != "unintelligible", is.analysis.group)
## Control:
## glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 1e+09))
##
## AIC BIC logLik deviance df.resid
## 1740.3 1842.0 -856.2 1712.3 10516
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.4174 -0.1399 -0.0842 -0.0538 19.3045
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subj (Intercept) 0.3870 0.6221
## trial.typeincon 2.0405 1.4285 0.78
## run 0.3482 0.5901 0.65 0.47
## trial.typeincon:run 0.8452 0.9193 -0.84 -0.94 -0.64
## Number of obs: 10530, groups: subj, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.2302 1.2276 -4.260 2.04e-05 ***
## trial.typeincon -0.1962 1.3758 -0.143 0.887
## run -0.3106 0.8119 -0.382 0.702
## trial.typeincon:run 0.9169 0.8861 1.035 0.301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trl.ty run
## tril.typncn -0.886
## run -0.938 0.858
## trl.typncn: 0.862 -0.953 -0.922
## convergence code: 0
## Model failed to converge with max|grad| = 0.403595 (tol = 0.001, component 1)
plot(rePCA(fit.error1.run)$subj)## get RT blups
blups <- as.data.frame(coef(fit1.het.trim))
blups %<>% rename(congr = "(Intercept)", stroop = "trial.typeincon") %>% tibble::rownames_to_column("subj")
## bind with error blups
blups %<>%
full_join(
data.frame(
subj = rownames(coef(fit.error1)$subj),
er.logit.stroop = coef(fit.error1)$subj$trial.typei, ## extract logits
er.logit.congr = coef(fit.error1)$subj[["(Intercept)"]]
) %>%
mutate(
er.logit.incon = er.logit.stroop + er.logit.congr, ## logit of error on incon trials
## blup stroop effect in units percent error:
stroop.er = (logit2prob(er.logit.incon) - logit2prob(er.logit.congr)) * 100
) %>%
dplyr::select(subj, stroop.er),
by = "subj"
)
## draw figure
plot.behav <-
blups %>%
mutate(subj = factor(subj, levels = subj[order(stroop, decreasing = TRUE)])) %>%
dplyr::select(subj, stroop, stroop.er) %>%
reshape2::melt() %>%
filter(!is.na(value)) %>%
ggplot(aes(subj, value)) +
facet_grid(
rows = vars(variable), scales = "free", switch = "y",
labeller = as_labeller(c(stroop = "Response time", stroop.er = "% error"))
) +
geom_segment(aes(x = subj, y = 0, xend = subj, yend = value), color = "grey50", size = geom.line.size) +
geom_point(fill = "grey30", color = "black", shape = 21, size = geom.point.size) +
coord_capped_cart(left = "both") +
xlab("Subject") +
theme(
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
strip.placement = "outside",
strip.background = element_blank(),
strip.text = element_text(size = axis.title.size),
axis.line.y = element_line(size = axis.line.size*0.6),
axis.text.y = element_text(size = axis.text.size),
axis.text.x = element_blank(),
axis.ticks.y = element_line(size = axis.line.size),
axis.ticks.x = element_blank(),
axis.title = element_text(size = axis.title.size),
axis.title.y = element_blank()
)## Using subj as id variables
ggsave(here("out", "behav", "stroop-blups.pdf"), height = 2.5, width = figwidth)
## write
write.csv(blups, here("out", "behav", "stroop_blups_rt_group201902.csv"), row.names = FALSE)
write.csv(behav, here("out", "behav", "behavior-and-events_group201902_with-subset-cols.csv"), row.names = FALSE)
behav.mod.objs <- list(
fit1.het.trim = fit1.het.trim,
er.stroopvar = er.stroopvar,
rt.stroopvar = rt.stroopvar,
rt.hom.v.het = rt.hom.v.het,
r.marginal.trim = r.marginal.trim
)
saveRDS(behav.mod.objs, here("out", "behav", "mod_objs.RDS"))behav %>%
mutate(subj = factor(as.factor(subj), levels = micinfo[order(micinfo$mic), "subj"])) %>%
ggplot(aes(subj, rt, color = mic)) +
geom_point(position = position_jitter(width = 0.1), size = 0.5, alpha = 0.3) +
geom_boxplot(position = position_nudge(1/3), notch = TRUE, width = 0.2) +
scale_color_brewer(type = "qual", palette = 2) +
theme(
axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0),
legend.position = c(0.8, 0.8)
)## Warning: Removed 111 rows containing non-finite values (stat_boxplot).
## Warning: Removed 327 rows containing missing values (geom_point).
behav %>%
group_by(mic, subj) %>%
summarize(freq = sum(rt == 0, na.rm = TRUE)) %>%
ggplot(aes(mic, freq)) +
geom_point(size = 2) +
labs(y = "frequency of rt == 0 per subj")## `summarise()` has grouped output by 'mic'. You can override using the `.groups` argument.
behav %>%
group_by(mic, subj, error) %>%
summarize(freq = sum(error, na.rm = TRUE)) %>%
ggplot(aes(mic, freq)) +
geom_point(size = 2, position = position_jitter(width = 0.1), alpha = 0.4) +
labs(y = "frequency of errors per subj")## `summarise()` has grouped output by 'mic', 'subj'. You can override using the `.groups` argument.
behav %>%
group_by(mic, subj, acc.final) %>%
summarize(freq = n()) %>%
ggplot(aes(mic, freq)) +
facet_grid(cols = vars(acc.final)) +
geom_point(size = 2, position = position_jitter(width = 0.1), alpha = 0.4) +
labs(y = "frequency of acc.final codings per subj")## `summarise()` has grouped output by 'mic', 'subj'. You can override using the `.groups` argument.
behav %>%
ggplot(aes(sample = rt, color = mic)) +
stat_qq(alpha = 0.3, size = 0.15) +
stat_qq_line(size = 0.25) +
facet_wrap(vars(subj)) +
scale_color_brewer(type = "qual", palette = 2)## Warning: Removed 111 rows containing non-finite values (stat_qq).
## Warning: Removed 111 rows containing non-finite values (stat_qq_line).
## Warning: Removed 216 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).
behav.mean.sd <- behav %>%
group_by(mic, subj, session) %>%
summarize(rt.mean = mean(rt, na.rm = TRUE), rt.sd = sd(rt, na.rm = TRUE))## `summarise()` has grouped output by 'mic', 'subj'. You can override using the `.groups` argument.
grid.arrange(
behav.mean.sd %>%
ggplot(aes(mic, rt.mean)) +
facet_grid(vars(session)) +
geom_point(position = position_jitter(width = 0.1)) +
geom_boxplot(position = position_nudge(1/3), notch = FALSE, width = 0.2),
behav.mean.sd %>%
ggplot(aes(mic, rt.sd)) +
facet_grid(vars(session)) +
geom_point(position = position_jitter(width = 0.1)) +
geom_boxplot(position = position_nudge(1/3), notch = FALSE, width = 0.2)
)## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing missing values (geom_point).
behav.rt.mic <- behav %>% filter(!is.implausible.rt, !is.artifactual.rt, acc == 1, mic != "unknown")
fit <- lmer(
rt ~ trial.type * mic + (trial.type | subj),
behav.rt.mic
)
summary(fit)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ trial.type * mic + (trial.type | subj)
## Data: behav.rt.mic
##
## REML criterion at convergence: 153769.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1277 -0.5433 -0.1395 0.3371 10.7427
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subj (Intercept) 13033 114.16
## trial.typeincon 1221 34.94 0.27
## Residual 27971 167.25
## Number of obs: 11737, groups: subj, 56
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 718.81 33.44 53.92 21.493 < 2e-16 ***
## trial.typeincon 91.14 12.28 52.63 7.425 9.73e-10 ***
## micmicro 106.70 37.73 53.94 2.828 0.00656 **
## trial.typeincon:micmicro -11.66 13.86 52.81 -0.841 0.40393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trl.ty micmcr
## tril.typncn 0.138
## micmicro -0.886 -0.122
## trl.typncn: -0.122 -0.886 0.137
fit.het <- lme(
rt ~ trial.type * mic,
random = ~ trial.type | subj,
weights = varIdent(form = ~ 1 | subj),
data = behav.rt.mic,
control = cl1
)
summary(fit.het)## Linear mixed-effects model fit by REML
## Data: behav.rt.mic
## AIC BIC logLik
## 150199.3 150663.6 -75036.66
##
## Random effects:
## Formula: ~trial.type | subj
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 114.64089 (Intr)
## trial.typeincon 32.60195 0.27
## Residual 158.98955
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | subj
## Parameter estimates:
## 115825 123117 130114 130518 132017 138837
## 1.0000000 1.3680098 0.9492074 0.7273461 1.1412919 0.9936276
## 141422 150423 158136 165032 171330 178243
## 1.0374676 1.3434331 0.6236905 1.6796577 0.9486992 1.5709797
## 178647 178950 197449 203418 204319 233326
## 1.1122407 0.7281896 1.2069358 0.7358493 1.6734081 2.0684527
## 300618 317332 346945 352738 393550 448347
## 0.5321995 0.6442487 0.8708450 1.6948459 0.7676649 1.1961969
## 580650 594156 601127 765864 814649 843151
## 2.0987110 0.6359147 1.1858526 0.5632233 0.9543389 1.2030043
## 849971 873968 877168 DMCC1328342 DMCC1596165 DMCC1971064
## 1.5060373 0.7877662 1.3172439 0.3994182 1.6274300 0.8102333
## DMCC2442951 DMCC2803654 DMCC2834766 DMCC3062542 DMCC4191255 DMCC5009144
## 0.4411166 0.6736317 0.6986990 0.9809815 0.6522699 0.5090851
## DMCC5195268 DMCC5775387 DMCC6371570 DMCC6418065 DMCC6661074 DMCC6671683
## 0.5882389 1.4022155 0.9880953 0.5389969 0.7096706 0.5231217
## DMCC6705371 DMCC6721369 DMCC7297690 DMCC8033964 DMCC8050964 DMCC9441378
## 0.7451601 0.5573345 0.7790957 0.4987241 0.8518175 1.3639635
## DMCC9478705 DMCC9953810
## 0.8554145 0.9115002
## Fixed effects: rt ~ trial.type * mic
## Value Std.Error DF t-value p-value
## (Intercept) 721.6388 33.45750 11679 21.568817 0.0000
## trial.typeincon 86.3618 11.09771 11679 7.781950 0.0000
## micmicro 102.5939 37.76622 54 2.716552 0.0088
## trial.typeincon:micmicro -5.8815 12.59782 11679 -0.466869 0.6406
## Correlation:
## (Intr) trl.ty micmcr
## trial.typeincon 0.165
## micmicro -0.886 -0.147
## trial.typeincon:micmicro -0.146 -0.881 0.162
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.8867575 -0.6220068 -0.1711832 0.4300865 10.0370354
##
## Number of Observations: 11737
## Number of Groups: 56
fit.het.ml <- update(fit.het, . ~ ., method = "ML")
fit.het.mic.ml <- update(fit.het.ml, . ~ ., weights = varIdent(form = ~ 1 | mic))
fit.hom.mic.ml <- update(fit.het.ml, . ~ ., weights = NULL)
(anova.mic <- anova(fit.het.ml, fit.het.mic.ml, fit.hom.mic.ml))## Model df AIC BIC logLik Test L.Ratio p-value
## fit.het.ml 1 63 150227.7 150692.0 -75050.84
## fit.het.mic.ml 2 9 153774.8 153841.2 -76878.42 1 vs 2 3655.162 <.0001
## fit.hom.mic.ml 3 8 153813.9 153872.8 -76898.93 2 vs 3 41.017 <.0001
modobjs <- list(
anova.mic = anova.mic,
mic.model.means = fit.het
)
saveRDS(modobjs, here("out", "behav", "mod_objs_mic.RDS"))## initial fit
behav.vset.rt <- behav %>% filter(acc == 1, !is.na(rt), rt < 3000, rt > 250, !is.analysis.group)
fit.vset <- lme(
rt ~ trial.type,
random = ~ trial.type | subj,
data = behav.vset.rt,
weights = varIdent(form = ~ 1 | subj),
control = lmeControl(maxIter = 1e5, msMaxIter = 1e5, niterEM = 1e5, msMaxEval = 1e5),
method = "REML"
)
## trim and re-fit
behav.vset.rt$resid.p <- resid(fit.vset, type = "p")
behav.vset.rt$is.far.out <- farout(behav.vset.rt$resid.p)
fit.vset.trim <- update(fit.vset, subset = !is.far.out)
## model error
fit.error.vset <- glmer(
error ~ trial.type + (trial.type | subj),
behav %>% filter(response.final != "unintelligible", !is.analysis.group) %>% mutate(error = 1 - acc),
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)## boundary (singular) fit: see ?isSingular
summary(fit.error.vset)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: error ~ trial.type + (trial.type | subj)
## Data:
## behav %>% filter(response.final != "unintelligible", !is.analysis.group) %>%
## mutate(error = 1 - acc)
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+09))
##
## AIC BIC logLik deviance df.resid
## 680.3 711.4 -335.2 670.3 3661
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.4709 -0.1620 -0.1031 -0.0076 9.6958
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subj (Intercept) 17.87 4.227
## trial.typeincon 10.12 3.181 -1.00
## Number of obs: 3666, groups: subj, 17
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.961 4.715 -2.537 0.0112 *
## trial.typeincon 7.983 4.657 1.714 0.0865 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## tril.typncn -0.998
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## extract predictions
blups.vset <- as.data.frame(coef(fit.vset.trim))
blups.vset %<>% rename(congr = "(Intercept)", stroop = "trial.typeincon") %>% tibble::rownames_to_column("subj")
## bind with error blups
blups.vset %<>%
full_join(
data.frame(
subj = rownames(coef(fit.error.vset)$subj),
er.logit.stroop = coef(fit.error.vset)$subj$trial.typei, ## extract logits
er.logit.congr = coef(fit.error.vset)$subj[["(Intercept)"]]
) %>%
mutate(
er.logit.incon = er.logit.stroop + er.logit.congr, ## logit of error on incon trials
## blup stroop effect in units percent error:
stroop.er = (logit2prob(er.logit.incon) - logit2prob(er.logit.congr)) * 100
) %>%
dplyr::select(subj, stroop.er),
by = "subj"
)
## plot
plot.behav.vset <-
blups.vset %>%
mutate(subj = factor(subj, levels = subj[order(stroop, decreasing = TRUE)])) %>%
dplyr::select(subj, stroop, stroop.er) %>%
reshape2::melt() %>%
filter(!is.na(value)) %>%
ggplot(aes(subj, value)) +
facet_grid(
rows = vars(variable), scales = "free", switch = "y",
labeller = as_labeller(c(stroop = "Response time", stroop.er = "% error"))
) +
geom_segment(aes(x = subj, y = 0, xend = subj, yend = value), color = "grey50", size = geom.line.size) +
geom_point(fill = "grey30", color = "black", shape = 21, size = geom.point.size) +
coord_capped_cart(left = "both") +
xlab("Subject") +
theme(
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
strip.placement = "outside",
strip.background = element_blank(),
strip.text = element_text(size = axis.title.size),
axis.line.y = element_line(size = axis.line.size*0.6),
axis.text.y = element_text(size = axis.text.size),
axis.text.x = element_blank(),
axis.ticks.y = element_line(size = axis.line.size),
axis.ticks.x = element_blank(),
axis.title = element_text(size = axis.title.size),
axis.title.y = element_blank()
)## Using subj as id variables
plot.behav.vset## write and save
ggsave(here("out", "behav", "stroop-blups-validation.pdf"), height = 2.5, width = figwidth/2)
write.csv(blups.vset, here("out", "behav", "stroop_blups_rt_group201902_validation.csv"), row.names = FALSE)
saveRDS(fit.vset.trim, here("out", "behav", "fit1-het-trim_group201902_validation_set.RDS"))## initial fit
behav.raset.rt <- behav %>% filter(acc == 1, !is.na(rt), rt < 3000, rt > 250, subj %in% subjs.analysis.red)
fit.raset <- lme(
rt ~ trial.type,
random = ~ trial.type | subj,
data = behav.raset.rt,
weights = varIdent(form = ~ 1 | subj),
control = lmeControl(maxIter = 1e5, msMaxIter = 1e5, niterEM = 1e5, msMaxEval = 1e5),
method = "REML"
)
## trim and re-fit
behav.raset.rt$resid.p <- resid(fit.raset, type = "p")
behav.raset.rt$is.far.out <- farout(behav.raset.rt$resid.p)
fit.raset.trim <- update(fit.raset, subset = !is.far.out)
## model error
fit.error.raset <- glmer(
error ~ trial.type + (trial.type | subj),
behav %>% filter(response.final != "unintelligible", !is.analysis.group) %>% mutate(error = 1 - acc),
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)## boundary (singular) fit: see ?isSingular
summary(fit.error.raset)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: error ~ trial.type + (trial.type | subj)
## Data:
## behav %>% filter(response.final != "unintelligible", !is.analysis.group) %>%
## mutate(error = 1 - acc)
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+09))
##
## AIC BIC logLik deviance df.resid
## 680.3 711.4 -335.2 670.3 3661
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.4709 -0.1620 -0.1031 -0.0076 9.6958
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subj (Intercept) 17.87 4.227
## trial.typeincon 10.12 3.181 -1.00
## Number of obs: 3666, groups: subj, 17
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.961 4.715 -2.537 0.0112 *
## trial.typeincon 7.983 4.657 1.714 0.0865 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## tril.typncn -0.998
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## extract predictions
blups.raset <- as.data.frame(coef(fit.raset.trim))
blups.raset %<>% rename(congr = "(Intercept)", stroop = "trial.typeincon") %>% tibble::rownames_to_column("subj")
## bind with error blups
blups.raset %<>%
full_join(
data.frame(
subj = rownames(coef(fit.error.raset)$subj),
er.logit.stroop = coef(fit.error.raset)$subj$trial.typei, ## extract logits
er.logit.congr = coef(fit.error.raset)$subj[["(Intercept)"]]
) %>%
mutate(
er.logit.incon = er.logit.stroop + er.logit.congr, ## logit of error on incon trials
## blup stroop effect in units percent error:
stroop.er = (logit2prob(er.logit.incon) - logit2prob(er.logit.congr)) * 100
) %>%
dplyr::select(subj, stroop.er),
by = "subj"
)
## plot
plot.behav.raset <-
blups.raset %>%
mutate(subj = factor(subj, levels = subj[order(stroop, decreasing = TRUE)])) %>%
dplyr::select(subj, stroop, stroop.er) %>%
reshape2::melt() %>%
filter(!is.na(value)) %>%
ggplot(aes(subj, value)) +
facet_grid(
rows = vars(variable), scales = "free", switch = "y",
labeller = as_labeller(c(stroop = "Response time", stroop.er = "% error"))
) +
geom_segment(aes(x = subj, y = 0, xend = subj, yend = value), color = "grey50", size = geom.line.size) +
geom_point(fill = "grey30", color = "black", shape = 21, size = geom.point.size) +
coord_capped_cart(left = "both") +
xlab("Subject") +
theme(
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
strip.placement = "outside",
strip.background = element_blank(),
strip.text = element_text(size = axis.title.size),
axis.line.y = element_line(size = axis.line.size*0.6),
axis.text.y = element_text(size = axis.text.size),
axis.text.x = element_blank(),
axis.ticks.y = element_line(size = axis.line.size),
axis.ticks.x = element_blank(),
axis.title = element_text(size = axis.title.size),
axis.title.y = element_blank()
)## Using subj as id variables
plot.behav.raset## write and save
ggsave(here("out", "behav", "stroop-blups-validation.pdf"), height = 2.5, width = figwidth/2)
write.csv(blups.raset, here("out", "behav", "stroop_blups_rt_group201902_analysis-reduced.csv"), row.names = FALSE)
saveRDS(fit.raset.trim, here("out", "behav", "fit1-het-trim_group201902_analysis-reduced_set.RDS"))